Roey Angel 2021-02-20
Here we attempt to detect ASVs that were labelled with 13C in our soil incubations using differential abundance modelling. Using DESeq2 (Love, Huber and Anders 2014) we compare the relative abundance of each ASV in the fractions where 13C-labelled RNA is expected to be found (>1.795 g ml-1; AKA ‘heavy’ fractions) to the fractions where unlabelled RNA is expected to be found (<1.795 g ml-1; AKA ‘light’ fractions). The method has been previously described in Angel et al., (2018).
set.seed(2021)
alpha_thresh <- 0.05
LFC_thresh <- 0.26
samples_prep_path <- "./"
data_path <- "./DADA2_pseudo/"
# Metadata_table <- "./AnCUE_Metadata_decontam.csv"
# Seq_table <- "DADA2.seqtab_nochim_decontam.tsv"
# Seq_file <- "DADA2.Seqs_decontam.fa"
Ps_file <- "Ps_obj_decontam_filt3.Rds"
Tree_file <- "./Tree/DADA2.Seqs_decontam_filtered.filtered.align.treefile"# Load phylogenetic tree
Tree <- read_tree(paste0(data_path, Tree_file))
# load and merge phyloseq object
readRDS(paste0(data_path, Ps_file)) %>%
merge_phyloseq(.,
phy_tree(Tree)
) -> Ps_obj_SIP
sample_data(Ps_obj_SIP)$Site %<>% recode_factor(Plesne = "Plešné",
Certovo = "Čertovo",
.ordered = TRUE)Because the DESeq2 models will be run on each gradient separately, we
need to subset This is easily done using
HTSSIP::phyloseq_subset (Youngblut, Barnett and Buckley
2018)
# split, ignore time points (for labelled ASV plots)
test_expr_1 <- "(Site == '${Site}' & Oxygen == '${Oxygen}' & Label..13C. == 'Unlabelled') | (Site == '${Site}' & Oxygen == '${Oxygen}' & Label..13C. == '${Label..13C.}')"
params_1 <- get_treatment_params(Ps_obj_SIP, c("Site",
"Oxygen",
"Glucose",
"Label..13C."),
"Label..13C. != 'Unlabelled'")
Ps_obj_SIP_noTime_l <- phyloseq_subset(Ps_obj_SIP, params_1, test_expr_1)
# names(Ps_obj_SIP_noTime_l) %<>%
# map(., ~str_remove_all(.x, "\\s\\|\\s.*")) %>%
# map(., ~str_remove_all(.x, "\\(|\\)|Site == |Hours == |Oxygen == |Label..13C. == |'")) %>%
# map(., ~str_replace_all(.x, "([0-9]+)", "\\1 h"))
# split, include time points (for DESeq2 modelling)
test_expr_2 <- "(Site == '${Site}' & Hours == '${Hours}' & Oxygen == '${Oxygen}' & Label..13C. == '${Label..13C.}') | (Site == '${Site}' & Hours == '${Hours}' & Oxygen == '${Oxygen}' & Label..13C. == '${Label..13C.}')"
params_2 <- get_treatment_params(Ps_obj_SIP, c("Site",
"Hours",
"Oxygen",
"Glucose",
"Label..13C."))
# Generate a list of subsetted phyloseq objects
Ps_obj_SIP_byTime_l <- phyloseq_subset(Ps_obj_SIP, params_2, test_expr_2)
names(Ps_obj_SIP_byTime_l) %<>%
map(., ~str_remove_all(.x, "\\s\\|\\s.*")) %>%
map(., ~str_remove_all(.x, "\\(|\\)|Site == |Hours == |Oxygen == |Label..13C. == |'")) %>%
map(., ~str_replace_all(.x, "([0-9]+)", "\\1 h")) Let us look first at the dissimilarity in community composition between the different fractions. If the labelling was strong enough we should see a deivation of (some of) the heavy fractions from the light ones. However, a lack of a significant deviation does not mean unsuccesful labelling because if only a small minority of the community was labelled we might not see it here (but we will, hopefully, see it using DESeq2 modelling).
(mod1 <- adonis2(vegdist(otu_table(Ps_obj_SIP), method = "horn") ~ Lib.size + Site * Oxygen * Hours ,
data = as(sample_data(Ps_obj_SIP), "data.frame"),
permutations = 999
))| Df | SumOfSqs | R2 | F | Pr(\>F) | |
|---|---|---|---|---|---|
| Model | 8 | 34.27409 | 0.6734282 | 75.0094 | 0.001 |
| Residual | 291 | 16.62085 | 0.3265718 | NA | NA |
| Total | 299 | 50.89494 | 1.0000000 | NA | NA |
plot_lib_dist(Ps_obj_SIP)Ps_obj_SIP %>%
scale_libraries(round = "round") ->
Ps_obj_SIP_scaled
plot_lib_dist(Ps_obj_SIP_scaled)(mod2 <- adonis2(vegdist(otu_table(Ps_obj_SIP_scaled), method = "horn") ~ Lib.size + Site * Oxygen * Hours,
data = as(sample_data(Ps_obj_SIP_scaled), "data.frame"),
permutations = 999
))| Df | SumOfSqs | R2 | F | Pr(\>F) | |
|---|---|---|---|---|---|
| Model | 8 | 34.61687 | 0.6748572 | 75.49891 | 0.001 |
| Residual | 291 | 16.67824 | 0.3251428 | NA | NA |
| Total | 299 | 51.29511 | 1.0000000 | NA | NA |
(mod3 <- adonis2(vegdist(otu_table(Ps_obj_SIP_scaled), method = "horn") ~ Lib.size + Site * Oxygen * Hours * Density.zone,
data = as(sample_data(Ps_obj_SIP_scaled), "data.frame"),
permutations = 999
))| Df | SumOfSqs | R2 | F | Pr(\>F) | |
|---|---|---|---|---|---|
| Model | 16 | 40.35622 | 0.786746 | 65.2535 | 0.001 |
| Residual | 283 | 10.93889 | 0.213254 | NA | NA |
| Total | 299 | 51.29511 | 1.000000 | NA | NA |
Site_disp <- betadisper(vegdist(otu_table(Ps_obj_SIP_scaled), method = "horn"), get_variable(Ps_obj_SIP_scaled, "Site"))
permutest(Site_disp)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.2375 0.237546 5.5678 999 0.016 *
## Residuals 298 12.7139 0.042664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(Site_disp)Oxygen_disp <- betadisper(vegdist(otu_table(Ps_obj_SIP_scaled), method = "horn"), get_variable(Ps_obj_SIP_scaled, "Oxygen"))
permutest(Oxygen_disp)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.03023 0.0302278 3.1681 999 0.087 .
## Residuals 298 2.84331 0.0095413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(Oxygen_disp)Hours_disp <- betadisper(vegdist(otu_table(Ps_obj_SIP_scaled), method = "horn"), get_variable(Ps_obj_SIP_scaled, "Hours"))
permutest(Hours_disp)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 4 0.2595 0.064877 4.7232 999 0.002 **
## Residuals 295 4.0520 0.013736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(Hours_disp)Density_disp <- betadisper(vegdist(otu_table(Ps_obj_SIP_scaled), method = "horn"), get_variable(Ps_obj_SIP_scaled, "Density.zone"))
permutest(Density_disp)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.34276 0.34276 33.288 999 0.001 ***
## Residuals 298 3.06845 0.01030
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(Density_disp)Ord <- ordinate(Ps_obj_SIP_scaled, "CAP", "horn",
formula = ~ Site * Oxygen * Hours * Density.zone)
explained <- as.numeric(format(round(eigenvals(Ord)/sum(eigenvals(Ord)) * 100, 1), nsmall = 1))
Ord_plt <- plot_ordination(Ps_obj_SIP, Ord, type = "samples", color = "Label..13C.", justDF = TRUE)
p_ord_joint <- ggplot(Ord_plt) +
geom_point(aes(
x = CAP1,
y = CAP2,
color = Label..13C.,
size = Density..g.ml.1.,
shape = Oxygen
), alpha = 2 / 3) +
guides(colour = guide_legend(title = "Labelling"),
size = guide_legend(title = "Density (g ml<sup>-1</sup>)"),
shape = guide_legend(title = "Oxygen")) +
scale_colour_locuszoom() +
# scale_colour_manual(values = Gradient.colours) +
# scale_fill_manual(values = Gradient.colours, guide = "none") +
labs(x = sprintf("CAP1 (%s%%)", explained[1]),
y = sprintf("CAP2 (%s%%)", explained[2])) +
coord_fixed(ratio = sqrt(explained[2] / explained[1])) +
theme(legend.justification = "top",
legend.title = element_markdown(size = 11)
) +
scale_size_continuous(breaks = round(c(seq(min(Ord_plt$Density..g.ml.1.),
max(Ord_plt$Density..g.ml.1.),
length.out = 5),
1), 4),
range = c(0.1, 5)) +
facet_grid(Site ~ Hours) +
# ggtitle("Joint analysis") +
NULL
save_figure(paste0(fig.path, "Oridnation"),
p_ord_joint,
pwidth = 10,
pheight = 8,
dpi = 600)
knitr::include_graphics(paste0(fig.path, "Oridnation", ".png"))Now run the differential abundance models using DESeq2. We then
filter the resutls to include only ASVs with Log_2_ fold change
>LFC_thresh and significant at
P<alpha_thresh. Lastly, we run ‘LFC-shrinking’ based on
Stephens (2016).
# generate a deseq object
DESeq_obj_SIP_byTime_l <- mclapply(Ps_obj_SIP_byTime_l,
function(x) {phyloseq_to_deseq2_safe(x,
test_condition = "Density.zone",
ref_level = "Light")},
mc.cores = nrow(params_2))
# run dds pipeline
DESeq_obj_SIP_byTime_l %<>% mclapply(.,
function(x) {DESeq(x,
test = "Wald",
fitType = "local")},
mc.cores = nrow(params_2)) # run dds pipeline
# extract results from a DESeq analysis
DESeq_res_SIP_byTime_l <- mclapply(DESeq_obj_SIP_byTime_l,
function(x) {
results(x,
altHypothesis = "greater",
alpha = alpha_thresh,
contrast = c("Density.zone", "Heavy", "Light"))}, # redundant if phyloseq_to_deseq2_safe() was used but doesn't hurt
mc.cores = nrow(params_2))
DESeq_res_SIP_byTime_LFC_l <- mclapply(DESeq_obj_SIP_byTime_l,
function(x) {
results(x,
lfcThreshold = LFC_thresh,
altHypothesis = "greater",
alpha = alpha_thresh,
contrast = c("Density.zone", "Heavy", "Light"))}, # redundant if phyloseq_to_deseq2_safe() was used but doesn't hurt
mc.cores = nrow(params_2)) # Extract results from a DESeq analysis
DESeq_res_SIP_byTime_LFC_shrink_l <- map(seq(length(DESeq_obj_SIP_byTime_l)),
~lfcShrink(DESeq_obj_SIP_byTime_l[[.x]],
res = DESeq_res_SIP_byTime_LFC_l[[.x]],
coef = "Density.zone_Heavy_vs_Light",
type = "ashr"))
names(DESeq_res_SIP_byTime_LFC_shrink_l) <- names(DESeq_res_SIP_byTime_LFC_l)
# Compare
plotMA(DESeq_res_SIP_byTime_l[[2]], ylim = c(-2,2))plotMA(DESeq_res_SIP_byTime_LFC_l[[2]], ylim = c(-2,2))plotMA(DESeq_res_SIP_byTime_LFC_shrink_l[[2]], ylim = c(-2,2))# summarise results (lfcShrink doesn't change the values)
# map2(DESeq_res_SIP_byTime_l, print(names(DESeq_res_SIP_byTime_l)), ~summary(.x)) # summarise results
# for (i in seq(1, length(DESeq_res_SIP_byTime_l))) { # didn't manage with map
# print(names(DESeq_res_SIP_byTime_l[i]))
# summary(DESeq_res_SIP_byTime_l[[i]])
# }
for (i in seq(1, length(DESeq_res_SIP_byTime_LFC_l))) { # didn't manage with map
print(names(DESeq_res_SIP_byTime_LFC_l[i]))
summary(DESeq_res_SIP_byTime_LFC_l[[i]])
}## [1] "Čertovo & 12 h & Anoxic & Labelled"
##
## out of 3026 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 0, 0%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 3, 0.099%
## low counts [2] : 7, 0.23%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 216 h & Anoxic & Labelled"
##
## out of 2878 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 37, 1.3%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 1484, 52%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 24 h & Anoxic & Labelled"
##
## out of 2755 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 1, 0.036%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 4, 0.15%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 48 h & Anoxic & Labelled"
##
## out of 2787 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 15, 0.54%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 1222, 44%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 216 h & Anoxic & Unlabelled"
##
## out of 2825 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 0, 0%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 5, 0.18%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 12 h & Oxic & Labelled"
##
## out of 3053 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 7, 0.23%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 8, 0.26%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 24 h & Oxic & Labelled"
##
## out of 2954 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 107, 3.6%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 905, 31%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 48 h & Oxic & Labelled"
##
## out of 2882 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 122, 4.2%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 331, 11%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 72 h & Oxic & Labelled"
##
## out of 2898 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 142, 4.9%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 1, 0.035%
## low counts [2] : 1329, 46%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Čertovo & 72 h & Oxic & Unlabelled"
##
## out of 2979 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 11, 0.37%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 6, 0.2%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 12 h & Anoxic & Labelled"
##
## out of 3031 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 27, 0.89%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 1678, 55%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 216 h & Anoxic & Labelled"
##
## out of 3210 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 0, 0%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 116, 3.6%
## low counts [2] : 3, 0.093%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 24 h & Anoxic & Labelled"
##
## out of 2848 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 22, 0.77%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 2340, 82%
## (mean count < 15)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 48 h & Anoxic & Labelled"
##
## out of 2758 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 12, 0.44%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 1527, 55%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 216 h & Anoxic & Unlabelled"
##
## out of 2862 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 0, 0%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 2, 0.07%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 12 h & Oxic & Labelled"
##
## out of 2923 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 57, 2%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 1, 0.034%
## low counts [2] : 1790, 61%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 24 h & Oxic & Labelled"
##
## out of 2548 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 101, 4%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 1061, 42%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 48 h & Oxic & Labelled"
##
## out of 2531 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 13, 0.51%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 1820, 72%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 72 h & Oxic & Labelled"
##
## out of 2556 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 115, 4.5%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 0, 0%
## low counts [2] : 1113, 44%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
##
## [1] "Plešné & 72 h & Oxic & Unlabelled"
##
## out of 2798 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.26 (up) : 3, 0.11%
## LFC < -0.26 (down) : 0, 0%
## outliers [1] : 7, 0.25%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
# Store labelled ASVs and save them to a file
# DESeq_res_SIP_byTime_l %>%
# map(., ~subset(.x, padj < alpha_thresh & log2FoldChange > LFC_thresh)) %>%
# map(., ~as.data.frame(.x)) %>%
# map(., ~rownames_to_column(.x, "ASV")) %>%
# bind_rows(., .id = "Comparison") %>%
# arrange(Comparison, desc(baseMean)) %T>%
# write_csv(., file = "DESeq2_byTime_a-0.05.txt") ->
# DESeq_res_SIP_byTime_df
# Store labelled ASVs and save them to a file
DESeq_res_SIP_byTime_LFC_shrink_l %>%
map(., ~subset(.x, padj < alpha_thresh & log2FoldChange > LFC_thresh)) %>%
map(., ~as.data.frame(.x)) %>%
map(., ~rownames_to_column(.x, "ASV")) %>%
bind_rows(., .id = "Comparison") %>%
arrange(Comparison, desc(baseMean)) %>%
separate(., "Comparison" ,c("Site","Hours", "Oxygen", "Label"), sep = " & ") ->
DESeq_res_SIP_byTime_LFC_sig_df
# grab the taxonomy of the ASVs
prune_taxa(DESeq_res_SIP_byTime_LFC_sig_df$ASV, Ps_obj_SIP) %>%
tax_table() %>%
as("data.frame") %>%
rownames_to_column("ASV") %>%
merge(., DESeq_res_SIP_byTime_LFC_sig_df, by = "ASV") %>% # watch out: this merge recycles values!
arrange(Site, Oxygen, Hours, baseMean) %>%
write_csv(., file = "DESeq2_byTime_a-0.05_LFC0-322.txt")DESeq_res_SIP_byTime_LFC_sig_df %>%
get_variable() %>%
select_if(is.numeric) %>%
vis_value()DESeq_res_SIP_byTime_LFC_sig_df %>%
get_variable() %>%
select_if(is.numeric) %>%
vis_cor()# DESeq_results <- DESeq_res_SIP_byTime_LFC_shrink_l[2]
# plot_DESeq(DESeq_results, Ps_obj_SIP, plot_title = names(DESeq_results))
DESeq_plots <- map(seq(length(DESeq_res_SIP_byTime_LFC_shrink_l)),
~plot_DESeq(DESeq_res_SIP_byTime_LFC_shrink_l[.x],
Ps_obj_SIP, plot_title = names(DESeq_res_SIP_byTime_LFC_shrink_l[.x])))
Certovo_DESeq <- ((DESeq_plots[[6]] +
theme(legend.position = "none") +
theme(axis.text.x = element_blank())) +
(DESeq_plots[[1]] +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank())) +
(DESeq_plots[[7]] +
theme(legend.position = "none",
axis.text.x = element_blank())) +
(DESeq_plots[[3]] +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank())) +
(DESeq_plots[[8]] +
theme(legend.position = "none") +
theme(legend.position = "none",
axis.text.x = element_blank())) +
(DESeq_plots[[4]] +
theme(legend.position = "none") +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank()) +
ylim(NA, 5)) +
(DESeq_plots[[9]] +
theme(legend.position = "none",
axis.text.x = element_blank())) +
(DESeq_plots[[2]] +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank())) +
(DESeq_plots[[10]] +
theme(legend.position = "none") +
ylim(NA, 5)) +
(DESeq_plots[[5]] +
theme(legend.position = "none",
axis.title.y = element_blank())) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = 'bottom'))
save_figure(paste0(fig.path, "Certovo_DESeq2"),
Certovo_DESeq,
pwidth = 8,
pheight = 10,
dpi = 600)
knitr::include_graphics(paste0(fig.path, "Certovo_DESeq2", ".png"))Plesne_DESeq <- ((DESeq_plots[[16]] +
theme(legend.position = "none") +
theme(axis.text.x = element_blank())) +
(DESeq_plots[[11]] +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank())) +
(DESeq_plots[[17]] +
theme(legend.position = "none",
axis.text.x = element_blank())) +
(DESeq_plots[[13]] +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank())) +
(DESeq_plots[[18]] +
theme(legend.position = "none") +
theme(legend.position = "none",
axis.text.x = element_blank())) +
(DESeq_plots[[14]] +
theme(legend.position = "none") +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank())) +
(DESeq_plots[[19]] +
theme(legend.position = "none",
axis.text.x = element_blank())) +
(DESeq_plots[[12]] +
theme(legend.position = "none",
axis.text.x = element_blank(),
axis.title.y = element_blank()) +
ylim(-3, NA)) +
(DESeq_plots[[20]] +
theme(legend.position = "none")) +
(DESeq_plots[[15]] +
theme(legend.position = "none",
axis.title.y = element_blank())) +
plot_layout(ncol = 2, guides = "collect") &
theme(legend.position = 'bottom'))
save_figure(paste0(fig.path, "Plesne_DESeq2"),
Plesne_DESeq,
pwidth = 8,
pheight = 10,
dpi = 600)
knitr::include_graphics(paste0(fig.path, "Plesne_DESeq2", ".png"))plot_combintions <- crossing(Site = c("Čertovo", "Plešné"),
Oxygen = c("Oxic", "Anoxic"))
Labelled_ASVs <- map(seq(length(Ps_obj_SIP_noTime_l)),
~plot_otus_by_density(Ps_obj_SIP_noTime_l[[.x]],
ASV2plot = filter(DESeq_res_SIP_byTime_LFC_sig_df,
Site == plot_combintions$Site[.x],
Oxygen == plot_combintions$Oxygen[.x])))
map(seq(length(Ps_obj_SIP_noTime_l)),
~save_figure(paste0(fig.path, "Labelled_ASVs_", paste(plot_combintions[.x, ], collapse = "_")),
Labelled_ASVs[[.x]],
pwidth = 16,
pheight = 12,
dpi = 600))## [[1]]
## [1] "05_Diff_abund_figures/Labelled_ASVs_Čertovo_Anoxic.svgz"
##
## [[2]]
## [1] "05_Diff_abund_figures/Labelled_ASVs_Čertovo_Oxic.svgz"
##
## [[3]]
## [1] "05_Diff_abund_figures/Labelled_ASVs_Plešné_Anoxic.svgz"
##
## [[4]]
## [1] "05_Diff_abund_figures/Labelled_ASVs_Plešné_Oxic.svgz"
plots2display <- list.files(path = paste0(fig.path),
pattern = "^Labelled_ASVs_(.*).png$",
full.names = TRUE)
knitr::include_graphics(plots2display)